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We study the charge conductivity of the one-dimensional repulsive Hubbard model at finite tem¬ 
perature using the method of dynamical quantum typicality, focusing at half filling. This numerical 
approach allows us to obtain current autocorrelation functions from systems with as many as 18 
sites, way beyond the range of standard exact diagonalization. Our data clearly suggest that the 
charge Drude weight vanishes with a power law as a function of system size. The low-frequency 
dependence of the conductivity is consistent with a finite dc value and thus with diffusion, despite 
large finite-size effects. Furthermore, we consider the mass-imbalanced Hubbard model for which 
the charge Drude weight decays exponentially with system size, as expected for a non-integrable 
model. We analyze the conductivity and diffusion constant as a function of the mass imbalance 
and we observe that the conductivity of the lighter component decreases exponentially fast with the 
mass-imbalance ratio. While in the extreme limit of immobile heavy particles, the Falicov-Kimball 
model, there is an effective Anderson-localization mechanism leading to a vanishing conductivity of 
the lighter species, we resolve finite conductivities for an inverse mass ratio of 77 > 0.25. 


I. INTRODUCTION 

The Hubbard model is a paradigmatic model in the 
theory of strongly correlated electrons, capturing some of 
the essential many-body effects due to short-range elec¬ 
tronic correlations in condensed matter physics: Mott- 
insulating behavior and the resulting localization of mag¬ 
netic moments with antiferromagnetic spin correlations. 
Moreover, the Hubbard model is the parent Hamiltonian 
for the Heisenberg and t-J model, which describe its low- 
energy physics in the strongly interacting regime [Tj-Q. 

The interest in the one-dimensional (ID) version of 
the model arises because of the existence of an exact so¬ 
lution based on the Bethc ansatz 3 and its relevance for 
quasi-ID materials [H-@, nanostructures IT0I IT2I and re¬ 
alizations with ultracold atomic gases in optical lattices 
d, [l3|] . A recent optical-lattice experiment has investi¬ 
gated the non-equilibrium c harg e transport in the two- 
dimensional Hubbard model 14]. 

The Hamiltonian of the ID repulsive Hubbard model 
is given by H = 1 hi with local terms 

hi = ( c la c l+i,a + h - c ‘) +U(ni ,t - - i) 

a 

(i) 

with CL+i,a = Ci.o-j where q. ct [c\ a ) annihilates (creates) 
a fermion with spin a =t, i on site l , and = cj CT c ; a is 
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the local density. L is the number of sites, fh is the hop¬ 
ping matrix clement, and U denotes the on-site Coulomb 
repulsion. 

Despite the success of the theory of such integrablc 
systems in computing many equilibrium properties, the 
quantitative and qualitative understanding of transport 
within linear response theory has proven to be a hard 
problem mm- While the zero-temperature transport 
properties are completely understood (see, e.g., [13), the 
main open questions concern transport of charge, spin, or 
energy at finite temperatures T > 0. The theory of the 
algebraic structure of the Bethe ansatz provides knowl¬ 
edge of local conservation laws, which can give rise to 
ballistic transport fl8i] . 

This ballistic transport is usually described via the 
Drude weight D , the zero-frequency contribution in the 
real part of the conductivity a(u>), 


Recr(w) = 2ttDS(uj) + <r reg (w). (2) 

As was argued by Zotos, Naef, and Prelovsek [Is|, a finite 
Drude weight exists if a lower bound is obtained from the 
Mazur inequality 
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where (•) is the thermodynamic average at temperature 
T. Such a bound exists if at least one conserved charge 
Qi has a finite overlap with the current operator j. The 
Qi = Yhi Qi,i are commonly ordered by their range, i = 1 
corresponding to particle number Q\ = N and i = 2 cor¬ 
responding to the Hamiltonian Q 2 = H. Q 3 has range 
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three (i.e., qi t 3 involves operators acting on three neigh¬ 
boring sites) and has the same structure as the energy- 
current operator, yet the two differ in the prefactor of 
one term fl 8 | . As a consequence, thermal transport in 
the one-dimensional Hubbard model is ballistic at any 
finite temperature T > 0 0 E3- Recently, it has been 
shown that there are also quasi-local conserved quantities 
in Bethe-ansatz integrable systems which can be crucial 
for some transport channels |20M22I] . Using the Mazur in¬ 
equality, one obtains a non-zero Drude weight for charge 
transport for any filling n = N/L ( N is the number of 
fermions) other than n = 1 / 2 , from considering only the 
leading non-trivial local conserved charge Q 3 of range 
three. The case of half filling has been discussed con¬ 
troversially, with some studies arguing in favor of a finite 
charge Drude weight D > 0 [Til l23l while others provided 
evidence for a vanishing D = 0 [2 ll T26I or at best a very 
small D [0 in the thermodynamic limit (tDMRG gives 
a small upper bound to the Drude weight). The situa¬ 
tion thus appears to be similar to spin transport in the 
spin-1/2 XXZ chain at zero magnetization, where also 
no local conservation law yields a non-zero bound to the 
spin Drude weight 0|, while numerical results [27l - l33j 
and Bethe-ansatz based calculations 0,0] strongly in¬ 
dicate a nonzero spin Drude weight at least in its gapless 
phase, with the possible exception of the point of full 
SU(2) symmetric exchange, i.e., the Heisenberg chain. 
For that model, though, quasi-local conservation laws 
have ultimately been identified as being at the heart of 
the ballistic spin transport mm at zero magnetization 
and in its gapless phase. 

The connection between (quasi-)local conservation 
laws and ballistic transport is closely related to how such 
conservation laws affect thermalization in integrable sys¬ 
tems [3jl]. Consider a quantum quench in which the force 
driving the current is turned off. If this initial condition 
leads to a finite value of (jQi), then the current will never 
completely decay back to zero. A simple example is the 
quench of a flux piercing a ring, which has been studied 
in this context [37( . 

Besides the question of the (divergent) zero-frequency 
contribution, the actual frequency dependence of the op¬ 
tical conductivity cy^fw) constitutes an equally inter¬ 
esting problem [26|, 1 38l l dOI . Some insight can be gained 
from effective low-energy theories such as bosonization 
[4ll443l ] . which is, however, limited to very low temper¬ 
atures and may not correctly capture effects due to in- 
tegrability without fine-tuning of parameters. An exact 
diagonalization study observed strong anomalous finite- 
size effects in cr reg (uj) of integrable Mott insulators IfiH . 
while many studies conclude that the dc conductivity 

a dc = lim cr reg (w) (4) 

uj —>-0 

is nonzero in such systems [2(1 HU- A recent density 
matrix renormalization group study suggests a generic 
divergence of a dc (T) at low temperatures with a dc oc 1 /T 
[2f| , different from the Fermi-liquid behavior a dc oc 1 /T 2 
that emerge in sufficiently high dimensions [44| . For the 


high-temperature regime, a lower bound for the diffusion 
constant V has been derived 0, reading 


T> > const. • 



( 5 ) 


(Note that T> ^ D.) While our primary interest is in 
the behavior in the linear response regime, we mention 
that numerical simulations of boundary-driven transport 
through open Hubbard chains also indicate diffusive high- 
temperature transport 0], 

In our work, we revisit the problem of charge trans¬ 
port in the Hubbard chain at half filling by employing 
the method of dynamical quantum typicality (DQT). Ba¬ 
sically, this approach uses single pure states that are 
constructed to yield typical thermal behavior at finite 
temperature to compute the time dependence of corre¬ 
lation functions. In the current context of transport, 
this method has recently been applied to the calcula¬ 
tion of the spin Drude weight in XXZ chains [331 and 
to transport in various non-integrable models 117 lOl. 
Since only a pure state needs to be propagated in the 
DQT method, any means of propagating the wave func¬ 
tion such as a forward integration or Krylov-space based 
approaches can be used, giving access to system sizes as 
large as L = 18, which is comparable to what can be 
reached for the ground state via Lanczos methods. 

We extract the Drude weight from the long-time be¬ 
havior of current autocorrelation functions and study its 
finite-size dependence. We observe a power-law decay 
with system size to zero, which we interpret in the frame¬ 
work of the eigenstate thermalization hypothesis applied 
to integrable systems 501 • Thus, our results confirm the 
predictions of Ref. (HLHa], he., a vanishing Drude weight 
D = 0 at finite temperatures. We further analyze the 
optical conductivity, for which our data suggest a finite 
cr dc . Depending on how the time-dependent data are con¬ 
verted to frequency, one either recovers the anomalous, 
system-size dependent fluctuations discussed in 0 or 
one obtains a smooth, diffusive-like low-frequency depen¬ 
dence. 

The Hubbard model can equivalently be formulated as 
a spin-1/2 model defined on a two-leg ladder: spin-up and 
spin-down fermions live on the two separate legs, where 
the exchange is of XY type along the legs, while on the 
rungs the Hubbard interaction translates into an Ising in¬ 
teraction. This reformulation is, on the one hand, useful 
for numerical implementations, and on the other hand, 
there are several natural ways of breaking the integrabil- 
ity that emerge in this picture. Transport in various spin 
Hamiltonians defined on spin ladders has in fact been 
intensely investigated [ 2 ^, [47|, 0, [5ll - l55j . 

Here, we consider the mass-imbalanced Hubbard model 
as an example of a non-integrable system. The local 
Hamiltonian now takes the form: 


hi = ~ J2 [*" ( 4 , 


a C l+l,a + h - C - 


+U{ni^~){ni A ~), 
( 6 ) 



3 


i.e., we introduce different hopping matrix elements t cT , 
a =t : I, for the two fermionic species. We define the 
inverse mass ratio as 


In the limit of rj = 0, also known as Falicov-Kimball 
model, one naturally obtains perfectly insulating behav¬ 
ior at any temperature due to an effective Anderson- 
localization mechanism. In this case, all the local den¬ 
sity operators of the heavy species become conserved 
quantities, i.e., [H,ni t i\ = 0- Thus, for a given random 
distribution of immobile spin-down fermions, via the in¬ 
teraction term Urii^ni.i, one effectively obtains a diag¬ 
onal disorder potential for the light fermions with local 
potentials ei = U nj j. drawn from a binary distribution 
ei = 0, U. The translational invariance of the original 
model at a given density of nj. = N±/L is restored by 
averaging over many random distributions of the heavy 
fermions. 

We are interested in the dependence of the conduc¬ 
tivities and oq(w) of the heavy and light species, 

respectively, as a function of the inverse mass ratio 77 . 
First, we compute the associated Drude weights, which 
vanish approximately exponentially fast with system size, 
as expected for a non-integrable model [H3, [H, [57j. For 
intermediate values of 77 , we observe a regular form of 
Cf(u>) and oq/w). The dc conductivity of the heavy com¬ 
ponent appears to simply vanish quadratically with t\, 
while the presence of the heavy fermions leads to an ap¬ 
proximately exponential decay of the dc conductivity of 
the light fermions as a function of decreasing 77 , which we 
are able to resolve for 77 > 0.25. 

The mass-imbalanced Hubbard model has recently at¬ 
tracted renewed interest in the context of many-body lo¬ 
calization [BM Us} since several authors have considered 
the possibility of man y-b ody localization in translation- 
ally invariant systems [60H62|. In our model, interactions 
could thus potentially lead to a non-trivial effect in the 
strongly mass-imbalanced regime. Recent work has sug¬ 
gested, though, that there likely is no mass-imbalance 
driven localization-delocalization transition in our model 
at a nonzero 77 , but a quasi many-body localized behavior 
with anomalous diffusion at small values of 77 f63i] . These 
results are based on exact diagonalization with L < 10. 
Our results suggest a finite, albeit exponentially small dc 
conductivity at least for 77 > 0.25. 

The plan of the paper is the following. Section |TT] sum¬ 
marizes the definitions and expressions of the conductiv¬ 
ity, the Drude weight, and current autocorrelation func¬ 
tions. In Section cm we provide a brief introduction to 
the DQT method and its application to the calculation 
of finite-temperature current autocorrelation functions. 
Section CY] contains our results for the integrable Hub¬ 
bard chain at half filling, while we present our data and 
the discussion of the mass-imbalanced model in Sec. El 
We conclude with a summary and an outlook in Sec. EH 


II. DEFINITIONS 

Using the Jordan-Wigner transformation, the mass- 
imbalanced Fermi-Hubbard model can equivalently be 
formulated as a spin- 1/2 model defined on a two-leg lad¬ 
der, 

hi= E -^^S^ + S^SI^ + US^JS) 

< 7 =T4 

where spin-up and spin-down fermions live on the two 
separate legs and the Hubbard interaction translates into 
an Ising interaction. Our numerical implementation is 
formulated in the spin language. 

We derive the charge current from the continuity equa¬ 
tion (Il|, leading to j =j^+j l and j a = i hi] in 

the Hubbard notation. In the spin notation, 

3° = 1,* - S?,*S? + 1 , J (9) 

1 

is the spin current in the first (a =t) or second (cr =|) 
leg. We correspondingly study the two current autocor¬ 
relation functions at inverse temperature (3 = 1 /T 

n n\ __ Re (ATW.v) _ ReTr{e~P H j a (t)j a } 

’ LZ LTr{e~P H } ’ (1 ° } 

where the time argument of ja(t) refers to the Heisen¬ 
berg picture, j a ( 0) = j a , and C a ( 0) = t„/2 in the high- 
temperature limit /3 —> 0. 

From the time dependence of C a (t) we determine the 
quantities 

C a {t 1 ,t 2 )= 1 ( d tC a {t) (11) 

*2 *1 Jt i 

in a time interval [ti,< 2 ] where C a (t ) has decayed to its 
long-time value C{t\ < t < t 2 ) « C(t 00 ) and is 
practically constant. Thus, the quantities C a (ti,t 2 ) ap¬ 
proximate the finite-size Drude weights of the two legs 
given by 

At = 77 - lim C a (0,t 2 ). (12) 

Z 7 T *2 —>00 

We determine the frequency-dependent optical conduc¬ 
tivity Re 0 \r,t max (w) via the finite-time Fourier transfor¬ 
mation 

1 p — /3 UJ max 

Rea CT , tma » =--- / d te** C a {t ). (13) 

w Jo 

Here, the choice of a particular t max implies a frequency 
resolution Soj ~ 7 r/< max . In the thermodynamic limit 
L —> 00 , Re o’ (Ti t max (a;) is a smooth function on an ar¬ 
bitrarily small scale Scu —> 0 and does not depend on 
the actual value of t max chosen, as long as it is large 
compared to the current relaxation time [48(. For any fi¬ 
nite L, however, it is important to find a reasonable t max 
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where finite-size effects are well controlled. In particular, 
for integrable systems, finding such a f max can be a subtle 
issue, as discussed later in detail. Note that, to leading 
order in /3, Recr cr (a;) oc (3 and that Recr cr (a;) = a a (uj) in 
the high-temperature limit. 

If we find a (t max , L ) region with no significant depen¬ 
dence on f max and L, we extract the dc conductivity a a ,dc 
as the low-frequency limit 

o- CTi dc = lim o- t (w). (14) 

In case of vanishing Drude weights, cr a ,dc/x is identical 
to the time-dependent diffusion constant 

jQ f'tm&x 

Pa(*max) = - d tC a (t) (15) 

X Jo 

with x being the static susceptibility and reading, at /3 —> 

0, 


X Tt{(E,-sy 2 }-(»fE,s,'„}) 2 1 

J =- L -= 4' (16> 

In the case of significant finite-size Drude weights, how¬ 
ever, P CT (t max ) may not depend on f max and L, while 
<7 <t, dc clearly does. Therefore, in such cases, the time- 
dependent diffusion constant provides a useful alterna¬ 
tive for extracting transport coefficients on the basis of 
finite systems. Beyond technical aspects, T>{t ) also has a 
clear physical interpretation: It directly yields informa¬ 
tion on how spatial variances of density profiles evolve in 
time f(Ml - [G7f for any finite L. 


III. DYNAMICAL QUANTUM TYPICALITY 
A. Concept 


the pure state \ip) is chosen according to the unitary in¬ 
variant Haar measure ffiBL ;69] and, according to typicality 
[toI fuil , a representative of the statistical ensemble. 

The pure state | ip), and each \4>n), correspond to the 
limit of high temperatures /? —> 0. We incorporate fi¬ 
nite temperatures /3 ^ 0 by introducing |= 
exp(— /3H/2) I’i/jn) and rewriting the current autocorrela¬ 
tion function in Eq. (fTUl) in the form [331 . I47L I68l . l69l . l76l . l77l ] 
(skipping the index a for clarity) 

rm = Rc 'L N ('>pN(P)\j{t) j\lpN(P)) 

+ e(|^}), (18) 

where e(|'*/')) is a statistical error resulting from the ran¬ 
dom choice of |t/>}. This error vanishes when sampling 
over several \iji) is performed, i.e., e = 0. 

However, the central advantage of Eq. m is not the 
vanishing mean error e = 0 but the knowledge about the 
standard deviation of errors E(e). This standard devia¬ 
tion is bounded from above by [33L [68l l69l . Iz3, 


E(e) < O 


f \/Rc (j(t)jj(t)j) \ 

y L^/dTs J 


(19) 


where c? e ff is the effective dimension of the Hilbert space. 
In the limit of high temperatures /? —> 0, d e g = 4 L is 
the full Hilbert-space dimension. Consequently, if the 
length L is increased, E(e) decreases exponentially fast 
with L. At arbitrary /3, d e ff = Tr{exp[— /3(H— Elo)]} is the 
partition function with ground-state energy Eq , reflecting 
the number of thermally occupied states, and also scales 
exponentially fast with L [32j. Therefore, while the error 
is exactly zero in the thermodynamic limit L —> oo, this 
error can be already very small at finite but large L and 
sampling is unnecessary, as is the case for all examples 
considered in our work. 


In this section we first introduce a very accurate ap¬ 
proximation of current autocorrelation functions. This 
approximation then provides the basis for the numerical 
technique used throughout our work. The central idea 
is to replace the trace operation Tr{*} = ]Tb(i| • |i) in 
Eq. (flOl) by a single scalar product {ip\ • \ip), where | ip) 
is a single pure state drawn at random. Since we aim at 
describing the current dynamics in the full Hilbert space, 
| ip) is drawn at random in the full basis. Conveniently, 
| ip) is randomly chosen in the eigenbasis of the particle 
number, 

d,N 

l'0)=X]l^ Jv ^ IV’jv) =^(a s + ?6 a )|s), (17) 

N s 

where s = s(N) is a label for the eigenstates with particle 
number N. The coefficients a s and b s are random real 
numbers. To be precise, these coefficients are chosen ac¬ 
cording to a Gaussian distribution with zero mean. Thus, 


B. Numerical implementation 

Most importantly, the approximation in Eq. m can 
be calculated without knowing the eigenstates and eigen¬ 
values of the Hamiltonian. This calculation is based on 
the two auxiliary pure states 

|<M/M)> = e-^ 4 -^ 2 |<M , (20) 

\vn(P, t)> = e~ lHt j e~ im/2 IV'iv) ■ (21) 

Both states are time- and temperature-dependent and 
the only difference between the two states is the addi¬ 
tional current operator j in the r.h.s. of Eq. m Using 
these states, the approximation in Eq. (El reads 

u TE w <<M/ 3 ,o)l<M£,o)} ' 1 j 

Apparently, the full time and temperature dependence 
in Eq. (1^1) results from the evolution of the pure states 
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only, i.e., there the current operator j is simply applied 
to the initial or time-evolved states. 

For, e.g., |$jv(/3, t)}, the 0 dependence is generated by 
an imaginary-time Schrodinger equation, 

Z dp) l<M/3 ’° )) = f (23) 

and the t dependence by the usual real-time Schrodinger 
equation, 

t^\* N (p,t))=H\$ N (fi,t)). (24) 

These differential equations can be solved by the use of 
straightforward iterative methods such as, e.g., Runge- 

Kutta 0, E3,!z3- We use a massively parallel imple¬ 
mentation of a Suzuki-Trotter product formula or Cheby- 
shev polynomial algorithm 0, l 79| , allowing us to study 
quantum systems with as many as 2 L = 36 lattice sites 
(L = 18 in the fermionic language), where the Hilbert- 
space dimension is d = (D(IC) 11 ). As compared to exact 
diagonalization, this dimension is larger by orders of mag¬ 
nitude. Yet, we do not exploit translation invariance of 
Hamiltonian and current. This symmetry adds momen¬ 
tum as a good quantum number and an additional layer 
of parallelization S3 • 

In practice, we use the Chebyshev polynomial algo¬ 
rithm to compute e~^ H ^ 2 \i/>n )■ The results of this algo¬ 
rithm are exact to at least ten digits. For the propaga¬ 
tion in real time, we mostly use a unitary, second-order 
product formula algorithm with a time step Stth = 0.02, 
which is sufficiently small to guarantee that the total en¬ 
ergy is conserved up to at least six digits. Occasion¬ 
ally, we have used the Chebyshev polynomial algorithm 
to compute the real-time evolution: no significant dif¬ 
ferences between these and the product-formula results 
were found. Most of the simulations were carried out on 
JUQUEEN, the IBM Blue Gene/Q located at the Julich 
Supercomputer Centre. A simulation of the largest sys¬ 
tem studied in the present paper (36 spins) required 3 
TB of memory, the computation was distributed over 
131,072 (MPI) processes, the total elapsed time to carry 
out 400 times steps was about 10 hours (1.3 million core 
hours). 


IV. RESULTS FOR THE HUBBARD MODEL 

This section contains our results for the charge trans¬ 
port in the ID Hubbard model, focusing at half filling. 
We consider infinite temperature 0 = 1/T 0 unless 

stated otherwise. First, we discuss the overall time de¬ 
pendence of the current autocorrelation function for var¬ 
ious values of 17/th- Second, we extract the Drude weight 
D from the long-time behavior of C(t). Finally, we dis¬ 
cuss the frequency dependence of the regular part and its 
zero-frequency limit. 



FIG. 1. (Color online) Real-time decay of the current auto¬ 
correlation function C(t) for (a) [7/th = 4, (b) [7/th = 8, (c) 
[7/th = 16 for various L = 9,11,13,15,18 and at infinite tem¬ 
perature 0th —> 0 (solid curves and circles). For the largest 
L = 15 and 18, convergence to system-size independent val¬ 
ues is reached at times tth ~ 5. For comparison, tDMRG 
data from [2f| are included in (a) and (b) (dashed curves). 
The inset in (a) shows the L = 9 result for two different ini¬ 
tial random states (solid curve: first state, squares: second 
state), which demonstrates small statistical errors for this L 
already. 


A. Time dependence of autocorrelation functions 

Figures [Ha)- (c) show typical results for the real-time 
decay of the current autocorrelation function C[t) for 
U/th = 4,8,16, respectively, and several system sizes 
L < 18. The figures show C(t) for times up to tth < 8, 
where the dominant decay of C(t) from its initial value 
occurs. Typically, the data from these different L co¬ 
incide for tth < 2.5. Beyond tth = 2.5, C(t) is a 
monotonically decreasing function of system size as in¬ 
dicated by the arrow in Fig. [Lib). The figures further 
include real-time density matrix renormalization group 
(tDMRG) data from [20] for comparison. Our DQT re¬ 
sults are in excellent agreement with the tDMRG data. 

As U/th increases, C(t) approaches small values in¬ 
creasingly faster as a function of time. On the other 
hand, the larger U/th , the more high-frequent and pro¬ 
nounced are the oscillations in C(t). These are inherited 
from the large U/th limit, in which the spectrum con¬ 
sists of bands of eigenstates separated by gaps of order 
U. These bands correspond to excitations with multi¬ 
ple doublons. Thus, the oscillatory dynamics in C(t) at 
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FIG. 2. (Color online) (a) Long-time limit of the current 
autocorrelation function C(t ) for different U/th = 4,8,16, 
fixed L = 15, and high temperatures f3th —> 0. (b) Finite-size 
scaling of the Drude weight C(ti,t 2 ), as extracted from the 
time interval [tith,f 2 fh] = [12.5,25], in a log-log plot. As a 
guide to the eyes, power laws (dashed lines) and a function 
oc 1/L (solid line) are indicated. The inset in (b) shows, for 
U/th = 16 and L = 16, that C{t\,t 2 ) does not depend on the 
specific choice of t 2 - 



FIG. 3. (Color online) The same information as shown in Fig. 
[2] (b) but for fixed U/th = 8 and different /3 th = 0,0.1, 0.5. 
Data for /3th = 0.0 and 0.1 almost coincide. 


large U/th is quite similar to the behavior in the spin-1/2 
XXZ chain in the strong Ising limit [§7} and spin-1/2 XX 
ladders in the strong rung-coupling limit [47|. 


B. Drude weight 

In order to extract the non-decaying portion of C(t), 
which equals the Drude weight, much longer times than 
tth ~ 12 need to be considered [33j. Therefore, we dis¬ 
play C(t) for tth < 25 in Fig. CHa) for the example of 
U/th = 4,8,16 and for L = 15. At times tth > 10, the 
oscillations in C(t) have decayed to a sufficiently small 
amplitude and hence we estimate the Drude weight by av¬ 
eraging C(t) in the time window t £ [fi th = 12.5, t 2 th = 
25], yielding C(ti,f 2 )- Note that C(ti,t 2 ) does not de¬ 
pend on this specific choice of t 2 , see the inset of Fig. 

n.b). 

The resulting, L-dependent C{t\,t 2 ) are shown in 
Fig. [5](b) in a log-log plot. The system-size dependence 
of C(ti,t 2 ) is consistent with a 1/L decay of the Drude 
weight to zero as system size increases. This scaling of D 
with system size is typical for integrable systems: it has 
been observed for the spin Drude weight of the spin-1/2 
XXZ chain as well [29l. l30l. l32l. l33l]. Moreover, the Drude 
weight approximately measures the fluctuations of diag¬ 
onal matrix elements of the associated current operator 
[50l ]. Such system-size dependent fluctuations are com¬ 
monly investigated to access the validity of the eigenstate 
thermalization hypothesis (80l482| . For integrable sys¬ 
tems, most numerical studies indicate a slow, power-law 
decay of these fluctuations [sol f57l. fS3| . Most notably, our 
data are consistent with a vanishing Drude weight D = 0 
at infinite temperature, in agreement with [25j]. 

In principle, if the infinite-temperature Drude weight 
vanishes, this does not necessarily imply that D(T ) = 0 
at any finite T. To see this, one can write the Drude 
weight in a high-temperature expansion 

D(T) = ^ + £l + ... (25) 

where D 1 is the infinite-temperature Drude weight stud¬ 
ied in Fig. [2](b). To substantiate that in the Hubbard 
model at half filling D(T) = 0 at any finite T, we have 
also computed D(T ) at T/th = 2,10, where D also seems 
to vanish as L increases. This is illustrated in Fig. [3] 


C. Optical conductivity 

Since the Drude weight appears to vanish as L —> 00 , 
all weight in Recr(o;) will ultimately be in the regular 
part <T reg (w). This optical conductivity has recently been 
studied using tDMRG 0], where a finite dc conductivity 
was observed that diverges as <7dc ~ 1/T as temperature 
decreases. 

We here first demonstrate that it is indeed possible 
to extract the dc conductivity from our time-dependent 
data for C(t). At infinite temperature, the dc conductiv¬ 
ity CTdc/x is simply equal to the integral V{t) over C(t) 
as defined in Eq. (usd, i.e., connected to the diffusion 
constant by an Einstein relation. 
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FIG. 4. (Color online) (a) Time-dependent diffusion constant 
T>(t) for U/th = 16, various L = 9,11,13,15,16, and high 
temperatures /3 th —► 0. A plateau is clearly visible at in¬ 
termediate times before the finite-size Drude weight yields a 
linear increase in the long-time limit. The plateau height is in¬ 
dependent of L and the plateau width increases with L. (This 
behavior is almost identical to the XXZ chain at A > 1.) (b), 
(c) Frequency dependence of the optical conductivity Rea (a;), 
as resulting from t max th = 4,100. (c) does not respect the 
“better” limit of L —> oo first and t max th —> oo afterward. 
Apparently, (c) shows strong finite-size effects at both, lo = 0 
and u ^ 0. However, in the thermodynamic limit L —> oo, 
(c) seems to approach (b). 


At large U/th, 2? increases quickly and then settles 
into a plateau, as is evident from the example presented 
in Fig. Oja). At large times, V further increases, which 
is due to both the non-zero Drude weight on finite sys¬ 
tems and other finite-size effects. Plotting data for T> for 
several system sizes clearly suggests that finite-size data 
gradually approach the plateau value at longer times as 
well, see Fig. SKa). 

The presence of such a plateau, following the reason¬ 
ing of [64j] , suggests a finite dc conductivity and diffusion 
constant. As shown in Fig. [5j the diffusion constant ex¬ 
hibits a peculiar behavior at T = oo: As U/th increases, 
it saturates at a U-independent value. This saturation 
results from the structure of the energy spectrum in the 
large- U/th limit: It consists of bands separated by U 
that have a band width given by th- Since we are tak¬ 
ing the limit U/th —> oo after taking the limit T —> oo, 


co. 


Q 



FIG. 5. (Color online) Time-dependent diffusion constant 
1 ?(t) for various U/th = 4,16,32,64, fixed L = 15, and high 
temperatures /3 th —>■ 0. Clearly, the plateau value of T>(t ) 
becomes independent of U in the limit of large U. 




CO/t, 

n 


FIG. 6. (Color online) The same information as shown in 
Fig- SI yet here for U/th = 4. Extracting Recr(u;) via the 
L-independent FT strategy is not applicable here, since T)(t) 
does not exhibit a clear plateau [see (a)], due to the long-time 
tail in C(t). This plateau will also not occur for L = 18. 
Thus, in (b), we present Rea(cu) obtained from a full-time 
FT, which thus results in strong finite-size effects at small 
frequencies. 


the dominant contribution to scattering comes from in¬ 
terband processes. This behavior appears to be generic 
for systems with an emergent ladder-like spectrum and 
has also been observed in the Ising regime of spin-1/2 
XXZ chains [67] and in spin-1/2 XX ladders [47| . The 
independence of the diffusion constant on U observed in 
Fig. [5] also unveils that the lower bound of (45|, as given 
in Eq. ©, is not exhaustive in the large-17 /th regime. 

For the purpose of computing Recr(o;), the existence 
of the plateau implies that the asymptotic behavior has 

















been reached. Moreover, the value of the plateau in T>(t) 
is independent of system size for the parameters of Fig. [I] 
Thus, we will compare two ways of computing Recr(o;): 
(i) the first version uses the full time dependence of C(t), 
up to and including times where we clearly observe finite- 
size effects (later dubbed full-time FT); (ii) In the second, 
we restrict the time window for the Fourier transforma¬ 
tion to times at which we have system-size independent 
data for C(f) (later referred to as L-independent FT). 

The results of both approaches are presented in 
Figs. Hb) and (c), respectively. The full-time FT re¬ 
solves the strong finite-size dependent structures that 
were known to exist from Ref. [3f|. The positions of 
these sharp peaks shift to smaller frequencies as system 
size increases. An extrapolation of <r reg (w) to zero fre¬ 
quency is thus difficult to control. 

The behavior of < 7 reg (u;) computed using the 21- 
independent FT strategy, by contrast, is a very smooth 
function that strongly resembles the optical conductivity 
of a typical diffusive system. This is clearly related to the 
fast initial decay of C(t) [see the data shown in Fig.QJb)], 
and the corresponding establishment of the plateau in 
the integrated quantity V(t), which consequently allows 
us to estimate the dc limit under the assumption that no 
additional time dependence emerges in C (t) at very long 
times and large systems. We thus propose that when¬ 
ever such a plateau is present in 1 7(f), the cleanest way 
of computing <r reg (uj) is the (.-independent FT, in line 
with the reasoning of Refs. [Hj, H3, H3|. 

Figure | 6 ] shows data for U/t^ = 4 as an example for a 
case, in which no clear plateau in T>(t ) can be resolved 
with the accessible system sizes. Here, we thus com¬ 
pute <r reg (w) from the full available time series of (7(f), 
which is shown in Fig. [61(b). The optical conductivity 
has a broad maximum at w/th ~ U /th and an additional 
low-frequency peak at w/fh ~ 1 whose position shifts to 
small frequencies as L increases. The data would suggest 
a small or vanishing dc conductivity, which we believe 
does not reflect the behavior of an infinitely large system 
[compare Fig.[4](b)], since the low-frequency finite-size ef¬ 
fects likely screen the correct low-frequency dependence. 


V. RESULTS FOR THE MASS-IMBALANCED 
CASE 


In this section, we present our results for the mass- 
imbalanced cases r] = t^/t^ < 1 , where the model is non- 
integrablc. We start with the case 77 = 0, the Falicov- 
Kimball limit, and discuss the emergence of Anderson 
localization in this limit. Then we turn to the case of 
77 ~ 1/2 and study both, Drude weight and optical con¬ 
ductivity. Finally, we summarize the scaling of the diffu¬ 
sion constant as a function of 77 in the rj region accessible 
to our numerical method. 



FIG. 7. (Color online) (a) Real-time decay of the current 
autocorrelation function Cf(t) of the light d-component for 
U/t -r = 4, strong imbalance 77 = t^/t-f = 0, L = 9,11,15, and 
high temperatures j3tf —>• 0. Since Cf(£) is highly oscillat¬ 
ing after the initial decay, also the time-dependent diffusion 
constant T>f(t) does so in (b). Consequently, the usual extrac¬ 
tion of a diffusion constant would depend on the specific point 
in time considered. However, the time average still yields a 
reasonable diffusion constant, (c) Finite-size scaling of the 
time average for U/t-f = 4,8, as resulting from the time in¬ 
terval [ti = [12.5, 75], in a semi-log plot. Apparently, 

the scaling is non-trivial, but the decrease is consistent with 
insulating behavior in the thermodynamic limit L —> 00 . 


A. Falicov-Kimball limit 

In the Falicov-Kimball limit 77 = 0 the model simplifies 
to 

hi =-t t ( 4 t c i+l t + h.c.)+t/(rq )t —-)(n u —i). (26) 

For this simplified model all commute with all local 
Hamiltonians hi and with each other, 

[hi,n k ,i\ = [ni : i, n k , 4 .] = 0 , (27) 

l,k = 1Each (rii^ — 1/2) is thus conserved and 
yields a good quantum number e; = ± 1 / 2 , with 2 L dif¬ 
ferent sequences 

e(ro) = (ei(m),... ,e L (m)) , (28) 

m = 1 ,..., 2 l . As a consequence, the full Hamiltonian 
H = Y^i hi can be rewritten as a sum of 2 L uncoupled 
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FIG. 8. (Color online) Frequency dependence of the optical 
conductivity Rea-t/w) for (a) f ma xtt = 25, (b) tmaxH = 50 
for U/t- 1 - = 4, strong imbalance r) = t^/tf = 0, different L = 
11,13,15, and high temperatures /3tf —> 0. (c) shows (a) for 
U/tf = 8. The overall structure is independent of t m ax and 
L. For the dependence of the dc limit Recr-|-(a> —¥ 0) on L and 
W, see Df(ti,t 2 ) discussed before. 


Hamiltonians H(m ) = hi(m ), where 

hi(m ) = -t t (4 it C; +l t + h.c.)+C/ ei{m)(ni A --) (29) 

and the U part becomes a site-dependent potential given 
by the sequence e(m). For many m, e(m) can be under¬ 
stood as a sequence of random numbers drawn from a 
binary distribution [—1/2,1/2]. Therefore, remarkably, 
many uncoupled Hamiltonians H(m) can be interpreted 
also as the single-particle, Anderson problem for on-site 
disorder of strength U. Note that translation invariance 
is typically broken for a given m but restored by sam¬ 
pling over m. Note further that all m contribute at finite 
temperatures. 

Due to the analogy to the single-particle, Anderson 
problem and the strict one-dimensionality of the lattice, 
one expects perfectly insulating behavior in the thermo¬ 
dynamic limit L —> oo at all temperatures. Early on, this 
expectation has been verified in numerical calculations of 
the optical conductivity [H, 86 ] for /?> 0 and values 
of U where the localization length does not exceed lattice 
sizes accessible. Yet, the high-temperature limit /3 f-j- —> 0 
has not been studied. 

In Fig.[7][a) we show our results for the time-dependent 
current autocorrelation function Cq-(i) for (3t-f 0 , 


U/t- f- = 4, and different L = 9,11,15. Clearly, Cf (t) de¬ 
cays rapidly on a rather short time scale tt^ ~ 1. After 
this initial decay Cj-(f) approaches zero from the negative 
side but still shows small oscillations. Note that these os¬ 
cillations are no finite-size effects since curves for L = 11 
and 15 are practically identical to each other for the long 
times tt-f ~ 75 depicted in the figure. This curve for 
Cq-(i) yields the time-dependent diffusion constant T>^(t) 
shown in Fig.[7](b). After the initial increase of T>^(t) we 
find a strong decrease related to the region where Cq(f) is 
negative. Necessarily, T>^(t) also shows small oscillations 
not related to finite-size effects, as evident from compar¬ 
ing L = 11 and 15 again. 

The long-time oscillations of T>^(t) indicate that the 
dynamical process cannot be described by a diffusion con¬ 
stant in the strict sense. However, to extract an effective 
diffusion constant, we average T>^(t) over the long-time 
interval [ti t-f, *2 if] = [12-5,75]. In Fig. [7]]c) we depict 
the resulting as a function of L for H/f-j- = 4, 8 

in a semi-log plot. Apparently, this time-averaged quan¬ 
tity decreases as system size increases and may eventu¬ 
ally become zero in the thermodynamic limit L —> oo. 
Note that the scaling for small L is partially related to 
tiny finite-size Drude weights D- j-, entering T>^(t) via the 
relation T>^[t) oc D^t in the long-time limit. 

Next we turn to the optical conductivity. Since C^(t) 
and T>^(t) do not become constant in the long-time limit, 
the finite-time Fourier transform necessarily depends on 
the specific time interval chosen. Thus, we show in 
Figs.[51[a) and (b) the Fourier transform of U/tf = 4 data 
for f max tf = 25 and 50, where times t < f max where con¬ 
sidered in the Fourier transformation. While Figs. [ 8 ](a) 
and (b) differ with respect to details, the overall struc¬ 
ture does not depend on the specific choice of t max . In 
particular, the limit ui —> 0 is consistent with a vanishing 
dc conductivity. Note that this limit coincides with T>^(t) 
evaluated at tt-j- = 25 and 50, respectively. Similarly, our 
results indicate a vanishing dc conductivity for U/t^ = 8 , 
as shown in Fig.(5]]c). The small negative spectral weight 
is an artifact of the finite-time Fourier transform used and 
depends on the specific choice of t max . 

To summarize, our /3 f-j- —> 0 results are consistent with 
the interpretation of the model in terms of the single¬ 
particle, Anderson problem in one spatial dimension. 


B. Intermediate imbalance 

Next we discuss the region 0 < rj < 1 , where the model 
still is non-integrablc but the interpretation of the model 
in terms of the single-particle, Anderson problem is not 
possible any more. In fact, in this 77 region, we deal with 
a many-particle problem. 

We start with intermediate imbalance 77 = 0.4. In 
Fig. [9](a) we depict our results for the time-dependent 
current autocorrelation function C a (t) for the light 
(tr ='[') and the heavy (a =J.) component for U/t^ = 4 
and L = 9,15, still in the high-temperature limit j3t^ —> 
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FIG. 9. (Color online) Real-time decay of the current auto¬ 
correlation function C CT (t) for (a) U/tf = 4, (b) U /tf = 8 for 
77 = ti/tf = 0.4, both components <7 =t,4-, two L = 9,15, 
and high temperatures /3t-f —» 0. (c) Finite-size scaling of the 
Drude weight C CT (ti,f 2 ), as extracted from the time interval 
[ti t 2 it] = [25,50], in a semi-log plot. As a guide to the 
eyes, exponentials (dashed lines) are indicated. 


FIG. 10. (Color online) Frequency dependence of the optical 
conductivity Re < 7 , 7 ( 01 ) for the (a) light component <7 =4, (b) 
heavy component <7 = 4 , for U/tf = 8, 77 = t^/tf = 0.4, and 
high temperatures /3t -7 —>■ 0, as resulting from different L = 
11, 13,16 and fmaxtf = 10, 20. The independence of L and 
l m ax is evident, (c) U dependence of Re< 7 CT (o>) for a large 
L = 15 and long tmaxt-f = 20. A peak at o i/t-f = U/t\ is 
clearly visible. 


0. In Fig. Etb) we additionally show results for U/t -[• = 8 . 
For both components, C a (t) decays fast on a time scale 
ff-|- ~ 1 but revivals appear afterward. While these re¬ 
vivals are equally pronounced for <7 =4 and only Cq(t) 
becomes negative in the time interval ttf ~ 2.5. How¬ 
ever, any revivals eventually disappear and C„ (t) decays 
fully to approximately zero for <7 =4 and 4-- When com¬ 
paring curves for L = 9 and 15, it is also evident that 
finite-size effects are small on the physically relevant time 
scale. Thus, we are able to obtain information on C a (t ) 
in the thermodynamic limit L —> 00 without invoking 
intricate extrapolations. 

It is also evident from Figs. [9](a) and (b) that Drude 
weights D a are small, i.e., there is no long-time satura¬ 
tion of C a (t) at a significant positive value. However, 
it is instructive to discuss the actual value of the Drude 
weights in more detail. In Fig.[9](c) we show the finite-size 
scaling of C a (ti,t 2 ), as extracted from the time interval 
[fi f- 7 , £2 if] = [25,50], for tr =f,| and U/t^ = 4,8 in a 
semi-log plot. Interestingly, C a is larger for <7 = 4 . and 
does not depend on U. In all cases, the finite-size scaling 
of C a is remarkably well described by a simple exponen¬ 
tial decrease over three orders of magnitude, with a rela¬ 
tive value C a /C < j( 0) < 10 ~ 3 at L = 15. This exponential 


decrease is expected for strongly non-integrable models 
[ 47 I [56j and, moreover, is in accord with the eigenstate 
thermalization hypothesis [5C)|, [57[ . 

Since finite-size effects are small and C a {t) decays to 
approximately zero, we can accurately determine the op¬ 
tical conductivity by Fourier transforming data for finite 
L and t. In Figs. El a) and (b) we show the finite-time op¬ 
tical conductivity Recr CT (w) at U/t -^ = 8 for the light and 
heavy component, respectively. As expected, Re < 7 , 7 ( 0 ;) 
does neither depend on t max nor L and is a smooth func¬ 
tion of frequency w. Similarly to the integrable case 
77 = 0 , we find a broad maximum at ui/t-f ~ U/t^ for 
both <7. In contrast, the position of the additional peak 
at low u> depends on <7 but is roughly independent of U, 
as shown in Fig. ante). Most importantly, the dc con¬ 
ductivity is finite and its actual value is, relative to the 
amplitude of the low-w peak, larger for the heavy compo¬ 
nent <7 = 4 .. As a function of 17, this dc conductivity de¬ 
creases but is still finite for all U depicted, see Fig. Me). 
Therefore, at 77 = 0.4, we can exclude the existence of an 
insulator in the high-temperature limit (3t^ —> 0 . 
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FIG. 11. (Color online) Time dependence of the diffusion con¬ 
stant T> a (t) for (a) a =f, (b) a =4 for various 77 = ti/t-f = 
0.7,..., 0.2, a single U/tf = 8 , fixed L = 14, and high temper¬ 
ature Ptf —> 0. Apparently, is very sensitive to varying 

77, in contrast to Dj.(t). For imbalance 77 < 0.6, a plateau 
of V a {t) can be already seen for the L depicted, (c) 77 scal¬ 
ing of the plateau value for both components, as extracted 
at the point ttf — 15, in a semi-log plot. As a guide to the 
eyes, an exponential (dashed line) is indicated. Note that 

TArX " &a,dc- 


C. Scaling of diffusion constant and dc conductivity 

We eventually discuss the scaling of transport coeffi¬ 
cients as a function of imbalance 77 = For the 77 

discussed below, extracting the dc conductivity <J a ,Ac as 
Re a a (w —> 0 ) for finite L is equivalent to determining 
the plateau value of the time-dependent diffusion con¬ 
stant T> a (t). Therefore, we focus on an analysis of T> a (t), 
which can be concisely summarized for various rj. 

In Fig. [TlTa) we show the time-dependent diffusion 
constant T>t(t) of the light component for different 77 = 
0 . 7 , ..., 0 . 2 , a single U/t^ = 8, and fixed system size 
L = 14 . In Fig. hum we show V±(t) of the heavy 
component for the same set of parameters. Several com¬ 
ments are in order. First, for both <7 =t, 4 -, a plateau 
of X> CT (t) is clearly visible at times tt- j- ~ 15 for imbal¬ 
ances 0.3 < rj < 0 . 6 . We have checked that the plateau 
values T> a x coincide with the dc conductivity cf. 

Fig. [TO] for 77 = 0 . 4 , even though not shown explicitly 
for all 77. Second, for 77 > 0 . 6 , T> a {t) oc D a t due to 
strong finite-size Drude weights D a in the vicinity of the 


integrable point 77 = 1 , cf. Fig. [6] These finite-size ef¬ 
fects prevent us from determining the diffusion constant 
in the thermodynamic limit L —> 00. Third, for 77 < 0 . 3 , 
2?-|-(f) of the light component develops the small oscilla¬ 
tions around zero discussed in the context of the Falicov- 
Kimball limit 77 = 0 . These oscillations prevent us from 
determining the diffusion constant with sufficiently high 
accuracy. Fourth, T>^(t) is much more sensitive to vary¬ 
ing 77 than Note, however, that we depict T>i(t)/t 2 

rather than T>^(t). In this way, we do not show the trivial 
scaling T> j_ (t) oc tf resulting from the static scaling of the 
current operator oc tj.. 

In Fig. fTTT c) we depict the 77 dependence of the plateau 
values visible for L = 14 , in a semi-log plot. While 
we find T>±/t 2 « const., we observe a decrease of V -f as 77 
decreases, consistent with a simple exponential function. 
If we assume that this scaling continues to small 77 beyond 
the 77 range accessible, this assumption would imply the 
absence of a diffusion-localization transition at 77 ^ 0, 
consistent with the conclusions of [ffij. However, based 
on our results in Fig. Illl c). we cannot exclude the onset 
of many-body localization and a sudden drop of T>^ to 
zero at finite but small 77, as suggested in previous works 
0,111. Nevertheless, we can constrain the existence of 
a possibly localized regime to 77 <C 0 . 25 . 


VI. SUMMARY AND OUTLOOK 

In this work we studied finite-temperature charge 
transport in the one-dimensional repulsive Hubbard 
model at half filling. Using the method of dynamical 
quantum typicality, we were able to access system sizes 
much larger than what can be reached with full exact 
diagonalization, and with no restriction on the accessible 
time scales. This allowed us to extract the finite-size de¬ 
pendent Drude weight from the time dependence of cur¬ 
rent autocorrelation functions. The analysis of the finite- 
size dependencies indicated a vanishing Drude weight in 
the thermodynamic limit, in agreement with [ 25 j . We 
further computed the optical conductivity and provided 
evidence that it is (i) a smooth function of w at low fre¬ 
quencies and in the thermodynamic limit and (ii) that 
the dc conductivity is indeed finite, the latter in agree¬ 
ment with [ 26 | . 

As an example of a non-integrable model, we consid¬ 
ered the mass-imbalanced Hubbard chain. This model 
has recently been discussed in the context of many- 
body localization in translationally invariant systems 
0 , 16 ll . [ 63 ]. We demonstrated the absence of a Drude 
weight for large L , as expected for a non-integrable sys¬ 
tem. Our results for inverse mass ratios of 77 > 0.25 
indicated a small dc conductivity, that appears to van¬ 
ish exponentially fast as a function of decreasing 77. At 
intermediate 77, the system is thus a normal diffusive con¬ 
ductor, while at small 77, the emergence of small long-time 
oscillations in the current autocorrelation function give 
rise to slightly anomalous transport, in line with the con- 
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elusions of Ref. [63j]. 

Extensions of our work comprise the study of 
finite-temperature charge and spin transport in one¬ 
dimensional strongly correlated electron systems. For 
instance, there is an intriguing prediction on the role of 
spin drag in one dimension, which has been claimed to 
give rise to diffusive spin transport, while charge trans¬ 
ports remains ballistic at finite temperature [871 ] . Such 
questions as well as other effects due to a coupling of the 
various transport channels in the Hubbard model and its 
variants constitute a rich playground for future work. 
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